# simulate a cohort

source("communitySim.R")
source("simCancer.R")
source("caseNumbers.R")
dyn.load("simCancer.so")

if(F) { # simulate some data.  if F, will just load it.
thecancers = c("Colon","Lung","Breast", "Prostate", "Stomach")

totalVar = 0.4899106
totalSD = sqrt(totalVar)
propGroupSD = 0.5

parameters = list(
	   geno=log(1), env=log(1), "geno:env"=log(1),
		 size=150000, probGeno=0.2,
		 probEnv=0.3, genoMiss= 0.1,
		 envMiss = 0.1, oneEnvPerCommunity =T,
     cancerMissRate=0.2, falseCancerRate=0.00045,
		 sdGroup=propGroupSD * totalSD,
		 sdPerson=sqrt(1-propGroupSD^2)*totalSD,
		 Enrollmentprobs = c(20000, 40000, 50000, 40000),
		 agerange = c(35, 69),
		 Followup = c(5, 10, 20, 30),
		 Ncommunity =  50,
		 EqualCommunity=T
)


sig = 0.001


populationData = getPopData()

deathFile = list(
  "F"=read.table("../data/MortalityNocancerFemale.txt", header=T),
  "M"=read.table("../data/MortalityNocancerMale.txt", header=T)
)
DiffCancer= list(
	"F"=lambdaDiffCancer(myfile="../data/IncidentFemaleDiffCancer.csv",
	 	agegroup),
   M=lambdaDiffCancer(myfile="../data/IncidentMaleDiffCancer.csv",
	 	agegroup)
	 )
# change simcohort and getDeath functions to take deathData and DiffCancer arguments
# instead of reading files, so that it's easier to use cvd death data
# if DiffCancer is missing (it will be for CVD), make all event types the same.


lambda =list(F=list(x=seq(30, 90, by=5), y= exp((-0.5)*0.585)*c(116.4,169.7,272.6,401.8,555.4,746.8,1003.2,1254.8,1528.7,1741.4,
1903.1,1950.2,1922.5) / 100000)  ,
M=list(x=seq(30,90, by=5), y=exp((-0.5)*0.585)*c(61.8,81.2,123.6,247.8,466.8,877.5,1403.7,2096.6,2598.4,2887.2,
3086.4, 3156.3, 2927.5)/ 100000) )


datamat = simCohort(parameters, lambda, populationData, deathFile=deathFile, DiffCancer)

save(datamat, file="datamat.RData")

} else {
load("datamat.Rdata")
}

Dcancer = "Colon"
Dfollowup= 20

thisCancer = rep(F, parameters$size)
thisCancer[grep(Dcancer, datamat$CancerType)] = T

	ageEndFollowup = datamat$Age + Dfollowup - datamat$Enrollment
		hadEvent = datamat$Event < ageEndFollowup
		datamat$thisEventType = hadEvent &  thisCancer

	 
		# censored failure time
		datamat$thisEventTime = datamat$Event
		datamat$thisEventTime[!hadEvent] = ageEndFollowup[!hadEvent]
	


datamatSmall = datamat[sample(dim(datamat)[1],10000),]

library(survival)
thesurv = Surv(datamatSmall$Age, datamatSmall$thisEventTime,datamatSmall$thisEventType)


library(INLA)
library(inlaDevel)

forInla = SurvToInla(thesurv)




theformula =  forInla ~ geno*env+  Gender + f(CensusDivision, model="iid",param=c(1,1))


# exponentials, this works
res = inla(theformula, data=datamatSmall, family="exponential")
res = inla(theformulaBig, data=datamat, family="exponential")


# weibull, this also works
res = inla(theformula, data=datamatSmall, family="weibull", 
control.fixed=list(prec=0.01, prec.intercept=0.1, mean=0),
  control.data=list(initial=1.1,param=c(.1, .1) ))


# weibull cure, 4 prior params, doesn't work
res = inla(theformula, data=datamatSmall, family="weibullcure", 
control.fixed=list(prec=0.01, prec.intercept=0.1, mean=0),
  control.data=list(initial=c(1.1,.1),param=c(.1, .1,.1,.1) ))



# coxph, this works.
thesurv = Surv(datamat$Age, datamat$thisEventTime,datamat$thisEventType)

		coxfitAfterYears=try(coxph(thesurv~geno*env+
			strata(Gender,na.group=TRUE) + frailty(CensusDivision), 
  		data=datamat, control=coxph.control(iter.max=20)) ) 
  		

